Graph Search Algorithms¶
Contents¶
Graph Search¶
A graph is a mathematical way to represent road networks. Graphs consist of 3 sets:
vertices/nodes
edges
a set representing relations between vertices and edges
The nodes represent intersections, and the edges represent the roads themselves. A route is a sequence of edges connecting the origin node to the destination node.
Searching is the systematic examination of states to find a path from the start state to the goal state. Mathematically speaking, a graph can be represented by \(G\), where
\(G=(V,E)\)
For a graph \(G\), vertices are represented by \(V\), and edges by \(E\).
Each edge is a tuple \((v,w)\), where
\(w\), \(v \in V\)
Weight can be added as a third component to the edge tuple.
Search algorithms can be broadly classified into deterministic algorithms and stochastic algorithms.
In the former, the search algorithm follows a rigorous procedure and its path and values of both design variables and the functions are repeatable. For the same starting point, the algorithm will follow the same path whether you run the program today or tomorrow.
In the latter, the algorithm always has some randomness and the solution is not exactly repeatable. Based on the availability of information about the search space (e.g. the distance from the current state to the goal), deterministic search algorithms can be broadly classified into blink/uninformed and informed search.

Graph Traversal Algorithms¶
Graph traversal refers to a process that traverses vertices of a graph following a certain order (starting from user-input sources). This category of graph search algorithms only seeks to find a path between two nodes, without optimizing for the length of the final route. In applications where the weight of edges in a graph are all equal (e.g. 1), BFS and DFS algorithms outperform shortest path algorithms like Dijkstra’s.
Breadth-first Search (BFS)¶
BFS is an algorithm where the traversal starts at a specified node (the source or starting node) and continues along the graph layerwise, thus exploring all exploring all of the the current node’s neighbouring nodes (those which are directly connected to the current node). If a result is not found, the algorithm proceeds to search the next-level neighbour nodes.
explored ← empty
found ← False
while frontier is not empty and found is False do
add node to explored
for child in node.expand() do
found ← True
Using BFS, search for the shortest path between The Equestrian Statue and the Bahen Centre. This example uses the same data as in Getting Started.
Note
This book uses the smart_mobility_utilities package for some operations, in order to simplify the process of visualizing graphs. You can find out more about downloading and installing the package here.
Let’s first find the largest connected component centered around our location, with a specified distance on each side. The reference point is the centre of the University of Toronto’s downtown campus.
import osmnx
reference = (43.661667, -79.395)
G = osmnx.graph_from_point(reference, dist=500, clean_periphery=True, simplify=True)
To plot the network, we will also need to highlight the starting and ending nodes. For the sake of simplicity, we will use the node id directly. For more information on finding the closest node to a given coordinate, refer back to the Getting Started section’s example.
highlighted = [1907446268, 1633421938]
# marking both the source and destination node
nc = ['red' if node in highlighted else '#336699' for node in G.nodes()]
ns = [50 if node in highlighted else 8 for node in G.nodes()]
fig, ax = osmnx.plot_graph(G, node_size=ns, node_color=nc, node_zorder=2)
Let’s visualize the above graph on a ipyleaflet map, using a helper function from the smart_mobility_utilities package.
from smart_mobility_utilities.viz import draw_map
draw_map(G,highlight=highlighted, force_leaflet=True)
Warning
For the purposes of this map, we use the force_leaflet option so that the map will be rendered by ipyleaflet. Normally, when there are more than 1,000 nodes in a graph, ipyleaflet performance is very slow. The visualization tools in smart_mobility_utilities will automatically switch to folium when there are more than 1,000 nodes, unless the force_leaflet flag is used. See the docs for smart_mobility_utilities for more information.
Currently, each node in the above graph is represented as a python dict with many attributes that are of no interest to us. This makes accessing certain properties of nodes overly complicated and verbose. To minimize this, we can use the Node class from smart_mobility_utilities.common to redefine the nodes, and only retain key information like parent, edge length from parent, and the node itself. You can view the structure of Node here.
from smart_mobility_utilities.common import Node, cost
from tqdm.notebook import tqdm
from collections import deque
from time import process_time
# First convert the source and destination nodes to Node
origin = Node(graph=G, osmid=1907446268)
destination = Node(graph=G, osmid=1633421938)
# Setup a progress bar using tqdm, max is V+E
time_complexity = len(G.nodes)+len(G.edges)
bar = tqdm(total=time_complexity)
# Process run time tracking for later analysis
start_bfs = process_time()
route = []
frontier = deque([origin])
explored = set()
found = False
while frontier and not found:
node = frontier.popleft()
explored.add(node)
for child in node.expand():
if child not in explored and child not in frontier:
if child == destination:
route = child.path()
found = True
frontier.append(child)
bar.update(1)
end_bfs = process_time() - start_bfs
bar.close()
print("Route: \n",route,"\n\n Cost:\n",cost(G,route))
Route:
[1907446268, 55808224, 55808416, 55808284, 1721866234, 389678268, 4953810915, 389678267, 24960090, 24960068, 1258698109, 389678145, 24960070, 24960073, 24960076, 24960080, 6028561924, 5098988924, 389678131, 6028562356, 854322047, 389677908, 24959560, 242413453, 749951161, 7311057931, 389678216, 389678215, 389678226, 1633421933, 1633421938]
Cost:
1385.116
Note
We used collections.deque instead of a list because of the time-complexity required for pop and append.
With pythonic list, these two operations take O(n) time, while collections.deque operates in O(1) time, as it is based on LinkedList.
Additionally, tqdm.notebook was used because we are rendering the progress in a Jupyter Notebook. For normal use, you can import directly from tqdm.
As you can see, the BFS algorithm managed to find the route from source to destination in under a second, having transversed 91% of the graph to find it. Let’s plot the route.
fig, ax = osmnx.plot_graph_route(G, route)
Let’s also plot it on a map.
from smart_mobility_utilities.viz import draw_route
draw_route(G, route)
Depth-first Search (DFS)¶
The DFS algorithm is a recursive algorithm that uses the idea of backtracking. It involves exhaustive searches of all the nodes by going as deep as possible into the graph. When it reaches the last layer with no result, it “backtracks” up a layer and continues the search.
explored ← empty
found ← False
while frontier is not empty and found is False do
add node to explored
for child in node.expand() do
found ← True
As you may have the noticed, the only difference between DFS and BFS is in the way that frontier works. Rather than working down layer by layer (FIFO), DFS drills down to the bottom-most layer and moves its way back to the starting node (LIFO).
Let’s implement this algorithm with our previous example.
time_complexity = len(G.nodes) + len(G.edges) # Maximum is V + E
bar = tqdm(total=time_complexity)
# Process run time for later analysis
start_dfs = process_time()
route = []
frontier = deque([origin])
explored = set()
found = False
while frontier and not found:
node = frontier.pop()
explored.add(node)
for child in node.expand():
if child not in explored and child not in frontier:
if child == destination:
route = child.path()
found = True
continue
frontier.append(child)
bar.update(1)
end_dfs = process_time() - start_dfs
bar.close()
print("Route: \n",route,"\n\n Cost:\n",cost(G,route))
Route:
[1907446268, 55808205, 55808194, 55808408, 55808414, 8711144452, 55808328, 55808437, 3210497979, 389677988, 1686556839, 389677984, 50885180, 36607322, 389677990, 389677993, 390545921, 60654129, 60654120, 50897854, 50897859, 389678001, 389678002, 2143434369, 390550470, 389678003, 390548860, 389678004, 771950946, 984911356, 728157228, 306721042, 389678005, 2143487625, 389678007, 5277943137, 2498969982, 389677902, 390545068, 390545043, 306725181, 390545044, 771931704, 775377001, 771950967, 8608123052, 771931728, 8608123055, 8608123068, 3554867351, 390545045, 390545047, 390545049, 390545050, 389678203, 390545078, 390545077, 24959544, 389678013, 389678205, 389678206, 389678207, 3996667046, 3996667045, 389678209, 389678210, 389678054, 389678175, 1432347915, 389678044, 389678043, 215726254, 3983181527, 389678211, 389678177, 24959555, 389678042, 389678184, 389678183, 389678216, 389678215, 389678214, 24959557, 389678218, 389678219, 773004741, 773004737, 5567060881, 5567060879, 1005007860, 1005007861, 24959565, 24959569, 389678241, 7311150229, 7311150221, 7311150222, 7311150218, 7311142107, 7311142109, 24959589, 390545033, 152659384, 389678238, 389678239, 389678240, 389678225, 389678245, 389678229, 729406374, 389678246, 29604723, 2143440970, 1458386384, 391188278, 389678247, 249989991, 391188296, 249990004, 389678248, 389678249, 6123553651, 6532307387, 6532307390, 969631968, 409731632, 4579468982, 969631975, 389678250, 394502545, 394502565, 394502563, 2809034239, 4380884143, 4380884142, 4376693531, 389678136, 389677925, 389678134, 2557539827, 389678133, 389677909, 749952029, 389677908, 389678222, 7311057936, 7311057937, 6028562355, 2557542523, 389677907, 239055729, 389678039, 389678040, 389677889, 389678220, 749951161, 393676412, 7311036242, 7311057933, 2557539817, 2557539816, 389678227, 1633422235, 1633421933, 1633421938]
Cost:
6092.076
draw_route(G,route)
It is very evident that the paths generated by our DFS and BFS implementations are not the most direct route. This is because both DFS and BFS are algorithms that can find routes between two nodes, but make no guarantees that they will return the shortest path. Additionally, DFS generally returns “deeper” results as it traverses the entire depth of the graph and works backwards to find a solution.
Shortest Path Algorithms¶
Dijkstra Search¶
Dijkstra’s algorithm allows us to find the shortest path between any two vertices of a graph. The algorithm creates a tree of shortest paths from the starting vertex, the source, to all other points in the graph. Dijkstra’s algorithm , published in 1959 and named after Dutch computer scientist Edsger Dijkstra, is the base of several other graph search algorithms commonly used to solve routing problems in popular navigation apps.

The following pseudocode and `python` implementation for Dijkstra's algorithm has been modified to work with our OSM data. This is because graphs generated from maps will naturally have self-loops and parallel edges.
Parallel edges may result in a route that is not the shortest available, as the route length depends heavily on which parallel edge was chosen when a particular path was generated. In the example below, the shortest path may be returned as 7, if the first edge connecting 0 and 1 is chosen when calculating that path.

Self-loops also cause trouble for the original Dijkstra algorithm. If a graph contains a self-loop, it may be the case that the shortest path to a node comes from itself. At that point, we would be unable to generate a route.

These two issues are generally easy but non-trivial to avoid. For parallel edges, we select the edge with the lowest weight (shortest length), and discard any other parallel edge.
With self-loops, we can ignore the loop entirely as negative-weight loops do not exist in routing problems (a road cannot have negative length), and positive-weight loops cannot be part of a shortest path.
unrelaxed_nodes ← empty
seen ← empty
// initialize the lists
for node in graph
add node to unrelaxed_nodes
shortest_dist[source] ← 0
found ← False
while unrelaxed_nodes is not empty and found is False do
add node to seen
if node is destination then
found ← True
distance ← shortest_dist[node] + length of edge to child
if distance < shortest_dist[child] then
child.parent ← node
This is what it looks like in python:
import math
time_complexity = int(len(G.nodes) + len(G.edges) * math.log(len(G.nodes)))
bar = tqdm(total = time_complexity) # Maximum time is E log V
# Process run time tracking for later analysis
start_dijkstra = process_time()
# Using a set here avoids the problem with self loops
seen = set()
shortest_dist = {osmid: math.inf for osmid in G.nodes()}
unrelaxed_nodes = [Node(graph = G, osmid = osmid) for osmid in G.nodes()]
shortest_dist[origin.osmid] = 0
found = False
while len(unrelaxed_nodes) > 0 and not found:
node = min(unrelaxed_nodes, key = lambda node : shortest_dist[node.osmid])
# relaxing the node, so this node's value in shortest_dist is the shortest distance between the origin and destination
unrelaxed_nodes.remove(node)
seen.add(node.osmid)
# if the destination node has been relaxed then that is the route we want
if node == destination:
route = node.path()
found = True
bar.update(1)
continue
# otherwise, let's relax edges of its neighbours
for child in node.expand():
# skip self-loops
bar.update(1)
if child.osmid in seen: continue
child_obj = next((node for node in unrelaxed_nodes if node.osmid == child.osmid), None)
child_obj.distance = child.distance
distance = shortest_dist[node.osmid] + child.distance
if distance < shortest_dist[child_obj.osmid]:
shortest_dist[child_obj.osmid] = distance
child_obj.parent = node
end_dijkstra = process_time() - start_dijkstra
bar.close()
print("Route: \n",route,"\n\n Cost:\n",cost(G,route))
Route:
[1907446268, 55808205, 55808194, 1907446267, 8699033082, 8699033084, 6542457312, 4953810914, 55808233, 299625330, 389677953, 7967019556, 7967019566, 4923076695, 55808571, 55808582, 389678112, 389678113, 389678146, 2143434862, 2143434860, 7311083158, 1258707987, 389678121, 50885147, 389678122, 389677906, 50885141, 389678180, 2143436415, 2143494216, 2143494214, 389678185, 1633421950, 389678184, 389678183, 389678216, 389678215, 389678226, 1633421933, 1633421938]
Cost:
1029.481
draw_route(G, route)
Hill Climbing¶
The idea of the algorithm is quite simple:
Starting with a known (non-optimized) solution to a function that is to be optimized, the algorithm checks the neighbours of that solution, and chooses the neighbour that is “more” optimized. The process is repeated until no “better” solution can be found, at which point the algorithm terminates.
While the algorithm works relatively well with convex problems, functions with multiple local maxima will often result in an answer that is not the global maximum. It also performs poorly when there are plateaus (a local set of solutions that are all similarly optimized).
neighbours ← children of current
while min(neighbours) < current do
neighbours ← children of current
Here, we introduce a few new ideas.
First, we treat the route between two nodes as a function, the value of which is the distance between the two nodes. Second, we generate “children” of this function.
The function¶
We need to define a function \(f\) that is our target for optimization.
\(f(x)\) gives us the length of a route for some given route \(x \in Y\), where \(Y\) is the set of all possible routes between two specific nodes.
How do we generate \(x\)? We could just generate random permutations between the two nodes, filtering for permutations that are feasible, and optimize \(f\) over these random, sparse permutations.
However, this method is not reproducible (because the permutations change every run).
Instead, we make a deterministic policy that generates a number of \(x \in Y\) by successively “failing” nodes between the source and destination nodes. We then find the shortest path between the nodes before and after the “failed” nodes.
By failing the nodes in a deterministic fashion, we can say that we have a function and neighbourhood with defined size for a certain value so we can “rigorously” conduct a local search.
To generate our initial route and children routes, we will use the smart_mobility_utilities package. You can see how these routes are generated by consulting the documentation for that package.
Note
The implementation below uses a function called get_children which is located in the smart_mobility_utilities library. It offers both normal and multiprocessed versions in order to speed up processing time.
We will be using the non-multiprocessed version for this example. To see an implementation of the multiprocessed version, see the analysis at the end of this section.
from smart_mobility_utilities.common import randomized_search
from smart_mobility_utilities.children import get_children
import matplotlib.pyplot as plt
# Visualize the costs over time
costs = []
# Set the number of children to generate
num_children = 20
# Process run time for later analysis
start_hill = process_time()
current = randomized_search(G, origin.osmid, destination.osmid)
costs.append(cost(G,current))
print("Initial cost:",costs[0])
neighbours = get_children(G,current,num_children=num_children)
shortest = min(neighbours , key = lambda route : cost(G, route))
print("Initial min(children):",cost(G,shortest))
while cost(G, shortest) < cost(G, current):
current = shortest
neighbours = get_children(G,current,num_children=num_children)
shortest = min(neighbours , key = lambda route : cost(G, route))
costs.append(cost(G,current))
print(f"Current cost:",costs[-1],"|","min(children):",cost(G, shortest))
route = current
end_hill = process_time() - start_hill
plt.xticks(range(len(costs)))
plt.title("Hill climbing: Shortest Path")
p = plt.plot(costs)
print("Final cost:",costs[-1])
Initial cost: 1315.842
Initial min(children): 1045.913
Current cost: 1045.913 | min(children): 1040.293
Current cost: 1040.293 | min(children): 1045.913
Final cost: 1040.293
draw_route(G,route)
While the implementation above is deterministic in nature, the initial route is still randomized. That means that it’s possible to get different results across runs.
Hill climbing will generally return some decent results as there are few local optimal points in the route function. However, with larger search spaces that will naturally have more local maxima and plateaus, it will get stuck fairly quickly.
Beam Search¶
While Hill Climbing maintains a single “best” state throughout the run, beam search keeps \(k\) states in memory. At each iterations, it generates the neighbours for each of the \(k\) states, and puts them into a pool with the \(k\) states from the original beam. It then selects the best \(k\) routes from the pool to become the new beam, and this process repeats. The algorithm terminates when the new beam is equal to the old beam. As it is a local search algorithm, it is also susceptible to being stuck at local maxima.
A beam search with \(k=\infty\) is the same as a BFS. Because there is the risk that a state that would lead to the optimal solution might get discarded, beam searches are considered to be incomplete (it may not terminate with the solution).
beam ← random k routes from source to destination
add beam to seen
pool ← children of routes in the beam with consideration of seen + beam
last_beam ← nil
while beam is not last_beam do
beam ← the best k routes from pool
add beam to seen
pool ← children of routes in the beam with consideration of seen + beam
Here is the code in python:
import heapq
num_neighbours = 10
k = 10
seen = set()
# Process run time for later analysis
start_beam = process_time()
beam = [randomized_search(G,origin.osmid,destination.osmid) for _ in range(k)]
# the seen routes must be converted to a tuple to be hashable to be stored in a set
for route in beam: seen.add(tuple(route))
pool = []
for route in beam:
children = get_children(G,route,num_children=num_neighbours)
for child in children:
if tuple(child) in seen: continue
else:
pool.append(child)
seen.add(tuple(child))
pool += beam
last_beam = None
print(f"Cost of initial beam: {[cost(G, route) for route in beam]}" )
while beam != last_beam:
last_beam = beam
beam = heapq.nsmallest(k, pool, key = lambda route: cost(G, route))
for route in beam: seen.add(tuple(route))
pool = []
for route in beam:
children = get_children(G,route,num_children=num_neighbours)
for child in children:
if tuple(child) in seen: continue
else: pool.append(child); seen.add(tuple(child))
pool += beam
print(f"The cost of the beam: {[cost(G, route) for route in beam]}")
route = min(beam, key = lambda route : cost(G, route))
end_beam = process_time() - start_beam
print("-" * 100)
print(f"Route: {route} \n\nCost: {cost(G, route)}")
Cost of initial beam: [1265.343, 1952.928, 1419.695, 1123.258, 1588.986, 1748.207, 1273.673, 1193.022, 1579.068, 1280.335]
The cost of the beam: [1075.923, 1092.343, 1106.204, 1122.814, 1123.258, 1127.188, 1128.148, 1136.327, 1193.022, 1225.857]
The cost of the beam: [1059.313, 1065.575, 1075.923, 1092.343, 1100.737, 1106.204, 1122.814, 1123.258, 1127.188, 1128.118]
The cost of the beam: [1059.313, 1065.575, 1075.923, 1082.629, 1084.127, 1091.86, 1092.343, 1100.737, 1106.204, 1111.064]
The cost of the beam: [1059.313, 1065.575, 1075.25, 1075.733, 1075.923, 1082.629, 1084.127, 1086.559, 1087.519, 1091.86]
The cost of the beam: [1056.739, 1059.313, 1065.575, 1075.25, 1075.733, 1075.923, 1082.629, 1084.127, 1086.559, 1087.489]
The cost of the beam: [1056.739, 1059.313, 1065.575, 1070.435, 1073.793, 1075.25, 1075.733, 1075.923, 1082.629, 1084.127]
The cost of the beam: [1056.739, 1059.313, 1065.575, 1070.435, 1073.793, 1075.25, 1075.733, 1075.923, 1077.723, 1078.683]
The cost of the beam: [1056.739, 1059.313, 1065.575, 1070.435, 1073.793, 1075.25, 1075.733, 1075.923, 1077.723, 1078.653]
The cost of the beam: [1056.739, 1059.313, 1061.599, 1065.575, 1070.435, 1073.793, 1075.25, 1075.733, 1075.923, 1077.723]
The cost of the beam: [1056.739, 1059.313, 1061.599, 1065.575, 1070.435, 1073.793, 1075.25, 1075.733, 1075.923, 1077.723]
----------------------------------------------------------------------------------------------------
Route: [1907446268, 55808205, 55808194, 1907446267, 8699033082, 8699033084, 6542457312, 4953810914, 55808233, 299625330, 389677953, 7967019556, 7967019566, 4923076695, 55808571, 55808582, 389678111, 389678112, 389678113, 24960070, 24960073, 7311083154, 50885160, 389678121, 50885147, 389678122, 389677906, 50885141, 24959549, 389677891, 389678187, 389678186, 389678177, 24959555, 389678212, 389678213, 389678214, 389678215, 389678226, 1633421933, 1633421938]
Cost: 1056.739
draw_route(G,route)
A* Search¶
A* (pronounced A-star) search is an informed search algorithm widely used in pathfinding and graph traversal.
A* works by “greedily” choosing which vertex to explore next, based on a function:
\(f(V) = h(V) + g(V)\), where \(h\) is a heuristic, and \(g\) is the cost accrued up to that point.
PQ ← min heap according to A* Heuristic
A*-SEARCH(source,destination) return a route
explored ← empty
found ← False
while frontier is not empty and found is False do
add node to explored
for child in node.expand() do
found ← True
The Heuristic¶
The driving force behind A* is the selection of a new vertex (or node) to explore based on the lowest heuristic value. This heuristic value is computed by the following formula:
let \(dist(x,y)\) be a function that calculates the straight line distance between two nodes \(x,y\),
and let \(O\) be the origin node, and \(D\) be the destination node,
\(h(V) = dist(V,O) + dist(V,D)\) for any given node \(V\)
As the sum of the distance to the origin and destination is minimized when \(V\) lies on a straight line from \(O\) to \(D\), this heuristic prioritizes nodes which are “closer” to the straight-line distance from origin to destination.
Note
The implentation of the A* heuristic in smart_mobility_utilities defaults to calculating distances as if the Earth were flat. For local searches, this yields the best results. If the size of the search area is larger, it is better to calculate distance using the haversine_distance, which takes into account the curvature of the Earth.
This can be done by setting the distance function like so:
astar_heuristic(G,origin.osmid,destination.osmid, measuring_dist = haversine_distance)
# Get the A* Heuristic for all the nodes in the graph
from smart_mobility_utilities.problem import astar_heuristic
# Process run time for later analysis
start_astar = process_time()
toOrigin, toDestination = astar_heuristic(G, origin.osmid, destination.osmid)
route = []
frontier = list()
frontier.append(origin)
explored = set()
found = False
while frontier and not found:
# choose a node based on its heuristic value
node = min(frontier, key = lambda node : toOrigin[node.osmid] + toDestination[node.osmid])
frontier.remove(node)
explored.add(node)
# expand its children
for child in node.expand():
if child not in explored and child not in frontier:
if child == destination:
route = child.path()
found = True
continue
frontier.append(child)
end_astar = process_time() - start_astar
bar.close()
print("The route is ",route)
print("Route cost is", cost(G, route))
The route is [1907446268, 55808224, 55808227, 8699043521, 1907446267, 8699033082, 55808301, 6542457312, 4953810914, 55808233, 299625330, 389677953, 7967019556, 7967019555, 390547782, 389678107, 389678111, 389678112, 389678113, 389678146, 2143434862, 2143434860, 7311083158, 1258707987, 389678121, 50885147, 389678122, 1258707990, 7153999189, 389678124, 6028562355, 2557542523, 389677907, 389678038, 389678039, 389678040, 2143436381, 7311057930, 7311057931, 389678216, 389678215, 389678226, 1633421933, 1633421938]
Route cost is 1137.69
draw_route(G,route)
Bi-Directional Search¶
The purpose of bi-directional searches are to run two simultaneous, non-parallel searches; one starts at the origin and the other at the destination, with the goal of meeting somewhere in between.
This approach is more efficient because of the time complexities involved.
For example, a BFS search with a constant branching factor \(b\) and depth \(d\) would have an overall time complexity of \(O(d^b)\). By running two BFS searches in opposite directions with only half the depth (\(d/2\)), the time complexity becomes instead \(O((d/2)^b)+O((d/2)^b)\), which is much lower than the original \(O(d^b)\).
frontier_f ← initialized with source
frontier_b ← initialized with destination
explored_f ← empty
explored_b ← empty
found ← False
collide ← False // if front overlaps with back
found ← False
altr_expand ← False// expansion direction
while frontier_f is not empty and frontier_b is not empty and not collide and not found do
add node to explored_f
for child in node.expand() do
if child is destination then
found ← True
collide ← True
altr_expand ← not altr_expand
add node to explored_b
for child in node.expand() do
if child is origin then
found ← True
collide ← True
altr_expand ← not altr_expand
Let’s implement A* as a bi-directional algorithm. Our criteria for expansion, both forwards and backwards, will be the A* heuristic.
# Process run time for later analysis
start_bidirectional = process_time()
# define destination and origin for the backwards expansion
destination_b = origin
origin_b = destination
# get A*
toOrigin_f, toDestination_f = astar_heuristic(G, origin.osmid, destination.osmid)
toOrigin_b, toDestination_b = astar_heuristic(G, origin_b.osmid, destination_b.osmid)
route = []
f_value = lambda node: toOrigin_f[node.osmid] + toDestination_f[node.osmid]
b_value = lambda node: toOrigin_b[node.osmid] + toDestination_b[node.osmid]
frontier_f = list()
frontier_b = list()
frontier_f.append(origin)
frontier_b.append(origin_b)
explored_f = list()
explored_b = list()
collide = False
found = False
altr_expand = False # to alternate between front and back
while frontier_f and frontier_b and not collide and not found:
if altr_expand:
# remove node_f from frontier_f to expand it
node = min(frontier_f, key = lambda node : f_value(node))
frontier_f.remove(node)
explored_f.append(node)
for child in node.expand():
if child in explored_f: continue
if child == destination:
route = child.path()
found = True
break
# checking for collusion with the target expansion
if child in explored_b:
overlapped = next((node for node in explored_b if node == child))
# we don't take the overlapped node twice
route = child.path()[:-1] + overlapped.path()[::-1]
collide = True
break
frontier_f.append(child)
altr_expand = False
else:
# remove node_b from frontier_b to expand it
node = min(frontier_b, key = lambda node : b_value(node))
frontier_b.remove(node)
explored_b.append(node)
for child in node.expand():
if child in explored_b: continue
if child == destination_b:
route = child.path()[::-1] # we reverse the list because we expand from the back
found = True
break
if child in explored_f:
overlapped = next((node for node in explored_f if node == child), None)
route = overlapped.path()[:-1] + child.path()[::-1]
collide = True
break
frontier_b.append(child)
altr_expand = True
end_bidirectional = process_time() - start_bidirectional
print("The route is \n\n",route)
print("Cost of the route:",cost(G,route))
The route is
[1907446268, 55808224, 55808227, 8699043521, 1907446267, 8699033082, 55808301, 6542457312, 4953810914, 55808233, 299625330, 389677953, 7967019556, 7967019555, 390547782, 389678107, 389678111, 389678112, 389678113, 389678146, 2143434862, 2143434860, 7311083158, 1258707987, 389678121, 50885147, 389678122, 1258707990, 7153999189, 389678124, 6028562355, 2557542523, 389677907, 389678038, 389678039, 389678040, 2143436381, 7311057930, 7311057931, 389678216, 389678215, 389678226, 1633421933, 1633421938]
Cost of the route: 1137.69
draw_route(G,route)
Hierarchical Approaches¶
When facing routing problems at larger scales, such as those involving entire continents or graphs with millions of nodes, it is simply implausible to use basic approaches like Dijkstra or DFS. Instead, routing algorithms “prune” the search space in order to simplify the routing problem. Additionally, routing services may choose to precompute certain routes and cache them on servers, so that response times to user queries are reasonable.
Hierarchical search algorithms prune the search space by generating admissible heuristics that abstract the search space. Read more about the general approach of hierarchical methods at Faster Optimal and Suboptimal Hierarchical Search.
In this section, we’ll give a brief overview of two hierarchical approaches that aim to solve the shortest path problem, and show how their heuristics are computed. There will also be a python implementation of the Contraction Hierarchies example.
Highway Hierarchies¶
In this algorithm, the hierarchy “level” of each road/arc in the graph is calculated. This distinguishes the type of road segment (i.e. residential, national roads, highways). This is further supplemented by relevant data such as maximum designated driving speed, as well as number of turns in the road. After the heuristics are generated for the graph, the data is passed through a modified search function (bi-directional Dijkstra, A*, etc) that considers the distance to the destination and the potential expansion node class.
For example, the algorithm will generally consider highways as viable expansion nodes when it is still relatively further away from the target, and then will start to include national roads, and finally residential streets as it nears the destination.
While this approach “makes sense”, there are some disadvantages. First, the algorithm largely overlooks what kind of roads humans “prefer” to drive on. That is to say, which a highway might make sense for a given route, the user may “prefer” to take local roads (i.e. driving to a friend’s house who lives nearby). Secondly, highway hierarchies do not take into account factors such as traffic, which fluctuates often and adds significant cost to an “optimal” route.
You can learn more about contraction hierarchies here.
Contraction Hierarchies¶
While the highway hierarchies algorithm may be useful for speeding up shortest path searches, it only considers three levels of road hierarchies. On the other hand, the contraction hierarchies algorithm introduced by Contraction Hierarchies: Faster and Simpler Hierarchical Routing in Road Networks has the same number of hierarchies as nodes in the graph, which is beneficial as increased number of hierarchies results in more pruning of the search space.
Contraction¶
Taking any node \(V\) from a graph \(G\), remove \(V\) as well as any connected edges.
Add any number of edges to \(G\) such that the shortest distance for any pair of neighbouring nodes remains the same, even after \(V\) has been removed. Below is an example:

Suppose we are contracting Node 1. By removing the node, we have now changed the shortest path for 0→2, 0→3, and 2→3. We can restore it by creating new edges. After the creation of new edges, we then reinsert any removed nodes and edges. We then update that node with its hierarchy level (which is just the order of contraction, from \(1 to n\)). See the below updated graph:

This process is then repeated for all the other nodes. In our example, no other nodes need to be contracted, as removing any one of the remaining nodes does not affect the rest of the graph in terms of shortest path. This process will form shortcuts, which allow us to search the graph much faster, as we can ignore certain nodes that have been “pruned”.
Contraction Order¶
Any order for node contraction will result in a successful algorithm, however some contraction ordering systems minimizes the number of new edges added to the graph, and thus the overall running time.
To utilize this, we can employ the idea of edge difference (ED). The ED of a node is \(S-E\), where \(S\) is the number of new edges added if that node were to be contracted, and \(E\) is the number of edges that would be deleted. Minimizing ED is equivalent to minimizing the number of new edges added to the graph, and thus improves processing time. We can calculate the ED of a graph with the below function:
def calculate_edge_difference(G, shortest_paths):
edge_difference = list()
seenBefore = list()
for i in G.nodes():
# used in edge difference calculations
edges_incident = len(G[i])
# we will be deleting the node entry
# from the original shortest paths
# dictionary so we need to save its state
# for later iterations
contracted_node_paths = shortest_paths[i]
del shortest_paths[i]
# excluding the node that we have just contracted
new_graph = [*G.nodes()]
new_graph.remove(i)
# let's compute the new shortest paths between
# the nodes of the graph without the contracted
# node so we can see the changes and add arcs
# to the graph accordingly but that is in
# the algorithm itself
new_shortest_paths = dict()
for source in new_graph:
new_shortest_paths[source] = dict()
for destination in new_graph:
# path the contracted node "i" to compute new shortest paths accordingly
new_shortest_paths[source][destination] = dijkstra_with_contraction(G, \
source, \
destination, \
contracted = i)
# the add arcs to keep the graph all pairs shortest paths invariant
shortcuts = 0
for source in new_shortest_paths:
# we get a copy from the original and the new shortest paths dictionary
SP_contracted = new_shortest_paths[source]
SP_original = shortest_paths[source]
for destination in SP_contracted:
# this is statement so we don't add 2 arcs
# for the same pair of nodes
if [source, destination] in seenBefore: continue
seenBefore.append(sorted((source,destination)))
# if there is a difference between the original SP and
# post-contraction SP -- just add new arc
if SP_contracted[destination] != SP_original[destination]:
shortcuts += 1
# let's leave the dictionary as we took it
# from the last iteration
shortest_paths[i] = contracted_node_paths
# this is the value of the contraction
# heuristic for that node
ED = shortcuts - edges_incident
edge_difference.append((i, ED))
return edge_difference
Example #1: Bi-directional Dijkstra with contraction¶
After contracting a graph, we can run a bi-directional Dijkstra to compute all the shortest paths in the graph. Following this, any queries for the shortest path between any two nodes in the graph can be solved using a cached result.
There are some restrictions for our modified Dijkstra’s:
The forward expansion only considers arcs \(u,v\), where \(level(u) > level(v)\). This forms the upward graph, where only nodes with a higher contraction order can be relaxed.
The backward expansion only considers arcs \(u,v\), where \(level(u) < level(v)\). This forms the downward graph, where only nodes with a lower contraction order can be relaxed.
These restrictions help “prune” the search space and speed up the processing.
At the end, multiple paths are returned from source to destination, of which the shortest is selected to be the solution.
Let’s try this with the sample graph below:
import networkx
import random
import math
G = networkx.Graph()
# Computing each contraction level requires a new graph, which is costly in terms of memory.
# To avoid this, we use a flag on every node to show its contraction state.
# Add 14 nodes
for i in range(1,15):
G.add_node(i, contracted=False)
# Define the edges
edges = [
(1,2,{'weight':1}),
(1,3,{'weight':4}),
(2,3,{'weight':5}),
(2,4,{'weight':2}),
(3,4,{'weight':2}),
(3,7,{'weight':2}),
(3,8,{'weight':1}),
(3,9,{'weight':1}),
(4,5,{'weight':5}),
(5,10,{'weight':7}),
(6,7,{'weight':4}),
(6,8,{'weight':3}),
(6,10,{'weight':3}),
(6,5,{'weight':3}),
(6,9,{'weight':1}),
(7,8,{'weight':6}),
(8,9,{'weight':3}),
(8,13,{'weight':5}),
(9,12,{'weight':1}),
(9,10,{'weight':3}),
(10,11,{'weight':4}),
(11,12,{'weight':3}),
(11,13,{'weight':4}),
(12,13,{'weight':2}),
(14,1,{'weight':3}),
(14,13,{'weight':2})
]
# Add the edges to the graph and visualize it.
G.add_edges_from([*edges])
networkx.draw(G,with_labels=True)
Below is a modified Dijkstra, with the condition that nodes are only considered if contracted=False. Note that all flags are set to False at the end of the algorithm, as graphs are passed by reference.
def dijkstra_with_contraction(G, source, destination, contracted = None):
networkx.set_node_attributes(G, {contracted: True}, 'contracted')
shortest_path = dict()
heap = list()
for i in G.nodes():
if not networkx.get_node_attributes(G, 'contracted')[i]:
shortest_path[i] = math.inf
heap.append(i)
shortest_path[source] = 0
while len(heap) > 0:
q = min(heap, key = lambda node : shortest_path[node])
if q == destination:
networkx.set_node_attributes(G, {contracted: False}, 'contracted')
return shortest_path[q]
heap.remove(q)
for v in G[q]:
# if the node is contracted we skip it
if not networkx.get_node_attributes(G, 'contracted')[v]:
distance = shortest_path[q] + G[q][v]['weight']
if distance < shortest_path[v]:
shortest_path[v] = distance
networkx.set_node_attributes(G, {contracted: False}, 'contracted')
return math.inf # if we can't reach the destination
With this complete, the shortest path for every node pair in the graph can now be computed. This will be done without any contraction. Then, the ED heuristic can be calculated for the entire graph.
shortest_paths = dict()
for i in G.nodes():
# dictionary for every node
shortest_paths[i] = dict()
for j in G.nodes():
# that will be filled with the shortest path
# between it and other nodes
shortest_paths[i][j] = dijkstra_with_contraction(G, i, j)
edge_difference = calculate_edge_difference(G,shortest_paths)
edge_difference.sort(key=lambda pair: pair[1])
edge_difference
[(8, -5),
(10, -4),
(5, -3),
(7, -3),
(11, -3),
(14, 0),
(2, 2),
(1, 3),
(6, 4),
(13, 5),
(4, 7),
(12, 18),
(3, 26),
(9, 33)]
While this is a good start to optimizing the node contraction order, it is by no means perfect. Notice that the ED values calculated above assume the node is the only node removed from the graph. Because we are successively removing every node in the graph, the ED list will potentially become inaccurate after even the first contraction.
Recall however that any arbitrary contraction order results in a successful algorithm. While there are ED heuristics that are able to update the ED list after each contraction, each of them come with their own costs and benefits.
For the purposes of this example, our ED list will be enough.
# to keep track of the edges added after the algorithm finishes
edges_before = [*G.edges()]
current_graph = [*G.nodes()]
# iterating over the tuples (node, level)
# from the sorted edge difference list
for node_ED in edge_difference:
node = node_ED[0]
# now we will contract the given node through all iterations
networkx.set_node_attributes(G, {node: True}, 'contracted')
# we have already contracted the node
# so there is no need
new_graph = current_graph
new_graph.remove(node)
current_shortest_paths = dict()
for source in new_graph:
current_shortest_paths[source] = dict()
for destination in new_graph:
current_shortest_paths[source][destination] = dijkstra_with_contraction(G, \
source, \
destination)
for source in current_shortest_paths:
SP_contracted = current_shortest_paths[source]
SP_original = shortest_paths[source]
for destination in SP_contracted:
if source == destination: continue
if SP_contracted[destination] != SP_original[destination]:
print("we have added edge between ", source, destination," after contracting", node)
# it seems like we add two edges instead of one, but this is simple graph
# so adding edge from a to b and then adding edge from b to a are the same
# doesn't add a thing, we didn't include condition for that because it will
# complicate the algorithm
G.add_edge(source, destination, weight=SP_original[destination])
current_graph = new_graph
we have added edge between 1 13 after contracting 14
we have added edge between 2 13 after contracting 14
we have added edge between 13 1 after contracting 14
we have added edge between 13 2 after contracting 14
we have added edge between 1 4 after contracting 2
we have added edge between 4 1 after contracting 2
# new edges after adding additional arcs
edges_after = [*G.edges()]
print("# edges before", len(edges_before))
print("# edges after", len(edges_after))
# edges before 26
# edges after 29
We can visualize these “shortcuts” on the graph.
added_edges = list(set(edges_after) - set(edges_before))
# let's color these edges and draw the graph again
colors = ['r' if edge in added_edges else 'k' for edge in G.edges()]
networkx.draw(G, with_labels=True, edge_color=colors)
We can reformat the ED list to become a hierarchy, for the purposes of our bi-directional Dijkstra:
hierarchical_order = dict()
for order, node in enumerate(edge_difference):
hierarchical_order[node[0]] = order
hierarchical_order
{8: 0,
10: 1,
5: 2,
7: 3,
11: 4,
14: 5,
2: 6,
1: 7,
6: 8,
13: 9,
4: 10,
12: 11,
3: 12,
9: 13}
Let’s now find the shortest path from \(8\) to \(12\).
source = 8
destination = 12
# Generating the upward graph
# initializing
SP_s = dict()
parent_s = dict()
unrelaxed_s = list()
for node in G.nodes():
SP_s[node] = math.inf
parent_s[node] = None
unrelaxed_s.append(node)
SP_s[source] = 0
# dijkstra
while unrelaxed_s:
node = min(unrelaxed_s, key = lambda node : SP_s[node])
unrelaxed_s.remove(node)
if SP_s[node] == math.inf: break
for child in G[node]:
# skip unqualified edges
if hierarchical_order[child] < hierarchical_order[node]: continue
distance = SP_s[node] + G[node][child]['weight']
if distance < SP_s[child]:
SP_s[child] = distance
parent_s[child] = node
# Generating the downward graph
# initializing
SP_t = dict()
parent_t = dict()
unrelaxed_t = list()
for node in G.nodes():
SP_t[node] = math.inf
parent_t[node] = None
unrelaxed_t.append(node)
SP_t[destination] = 0
# dijkstra
while unrelaxed_t:
node = min(unrelaxed_t, key = lambda node : SP_t[node])
unrelaxed_t.remove(node)
if SP_t[node] == math.inf: break
for child in G[node]:
# skip unqualified edges
if hierarchical_order[child] < hierarchical_order[node]: continue
distance = SP_t[node] + G[node][child]['weight']
if distance < SP_t[child]:
SP_t[child] = distance
parent_t[child] = node
With these, we now merge the common settled nodes from SP_d and SP_s, and find the minimum sum of values, This is our shortest path: the collission between the forward and backward expansions.
minimum = math.inf
merge_node = None
for i in SP_s:
if SP_t[i] == math.inf: continue
if SP_t[i] + SP_s[i] < minimum:
minimum = SP_t[i] + SP_s[i]
merge_node = i
print("Minimum:",minimum)
print("Merge node:",merge_node)
Minimum: 3
Merge node: 9
# Trace the route back
def route_dijkstra(parent, node):
route = []
while node != None:
route.append(node)
node = parent[node]
return route[::-1]
route_from_target = route_dijkstra(parent_t, merge_node)
route_from_source = route_dijkstra(parent_s, merge_node)
route = route_from_source + route_from_target[::-1][1:]
route
[8, 3, 9, 12]
Let’s compare this to networkx’s built in solver:
from networkx.algorithms.shortest_paths.weighted import single_source_dijkstra
single_source_dijkstra(G,source, destination)
(3, [8, 3, 9, 12])
Pruning¶
We can actually check to see how many nodes were “pruned” from our graph, using the upward and downward graphs generated earlier.
unvisited = 0
for s_node, s_dist in SP_s.items():
for t_node, t_dist in SP_t.items():
if s_node == t_node and s_dist == t_dist == math.inf:
unvisited += 1
print(f"""Skipped {unvisited} nodes from a graph with {len(G)} total nodes,
resulting in pruning {unvisited/len(G)*100}% of the nodes in our search space.""")
Skipped 7 nodes from a graph with 14 total nodes,
resulting in pruning 50.0% of the nodes in our search space.
Summary¶
Let’s use the processing run times we collected earlier to show a comparison of the processing times for the different algorithms discussed here.
print("BFS:",end_bfs,"s")
print("DFS:",end_dfs,"s")
print("Djikstra:",end_dijkstra,"s")
print("Hill Climb:",end_hill,"s")
print("Beam Search:",end_beam,"s")
print("A* Search:",end_astar,"s")
print("Bi-directional A*:",end_bidirectional,"s")
BFS: 0.013275471000000039 s
DFS: 0.02681750300000374 s
Djikstra: 0.1345126840000006 s
Hill Climb: 15.583480621 s
Beam Search: 135.696005033 s
A* Search: 0.007053232000004073 s
Bi-directional A*: 0.013135845000022073 s